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ABSTRACT 

A method is described for preprocessing projection data prior to image reconstruction in single-photon emission 
computed tomography. The projection data of the desired spatial distribution of emission activity is blurred by the point- 
response function of the collimator that is used to define the range of directions of gamma-ray photons reaching the detector. 
The point-response function of the collimator is not stationary, but depends on the distance from the collimator to the point. 
Conventional methods for deblurring collimator projection data are based on approximating the actual depth-dependent point- 
response function by a spatially-invariant blurring function, so that deconvolution methods can be applied independently to 
the data at each angle of view. The method described in this paper is based on Fourier analysis of the multi-angular data set as 
a whole, using special depth-dependent characteristics of the Fourier coefficients to achieve spatially-variant inverse filtering 
of the data in all views simultaneously. Preliminary results are presented for simulated data with a simple collimator model. 

1. INTRODUCTION 

Nuclear medicine involves the determination of functional and metabolic information about selected human organ 
systems from the distribution of radioactive tracers administered to the patient The conventional instrument for imaging 
the distribution of activity combines a gamma camera to detect and map the received radiation and a collimator to limit the 
range of directions of the detected photons. A stationary gamma camera with its associated collimator produces a projection 
image of the activity distribution within its field of view. This type of image has many uses in nuclear medicine, but it 
has the disadvantage that overlying and underlying parts of the activity distribution appear superimposed, analogous to the 
superposition problem inherent in the conventional x-ray projection image. Single-photon emission computed tomography 
(SPECT) 1 is a method for eliminating the superposition problem in nuclear medicine imaging. SPECT produces images of 
cross-sectional slices free from out-of-slice contributions by computer processing of the set of projection data that is acquired 
by a collimated gamma camera for many angles of view as it rotates around the patient. 

Although SPECT has now found widespread use in clinical nuclear medicine 2 , there remain some significant technical 
challenges in the design of the instrumentation and in the processing of the data. The design of the collimator involves 
an inherent compromise between resolution and sensitivity. As the collimator holes are made narrower and/or longer, the 
resolution is improved, but fewer of the emitted photons are accepted by the collimator. Also, the acceptance angle of a 
collimator hole relative to a point source depends both on the distance from the point to the collimator and on the distance 
from the point to the line formed by extending the central axis of the collimator hole into the activity distribution . The 
point-response function of the collimator and the consequent blurring in the detected projection of a point source are therefore 
depth-dependent. As the collimator is rotated around the patient, a fixed point source that is offset from the center of rotation 
will contribute greater or lesser amounts to the projection data as the distance from the point to the collimator varies. Our 
approach to the problem of depth-dependent collimator blurring is the main subject of this paper. 

In addition to the problem of position-variance of the collimator's point-response function, there are other significant 
problems that limit the quantitative capabilities of SPECT 4 , Because the radiation dose to the patient must be minimized, 
and because the collimator accepts only a very small fraction of the emitted photons, the number of photons contributing 
to each measurement is relatively small, resulting in large statistical uncertainties in the measured data. Other significant 
problems that affect the quality of the data are caused by photon attenuation and photon scattering between the original 
point of emission and the detector. 

In this paper we describe the basic principles of our approach to the correction of depth-dependent collimator blurring 
in SPECT. We use a simplified model of the actual image-formation process in order to concentrate on the basic ideas. We 
consider a two-dimensional (2-D) analog of the actual 3-D situation: specifically, we consider a 2-D distribution of activity 
in the plane, with 1-D projection data collected by a 1-D collimator for each angle of view around the activity distribution. 
The present formulation of our approach does not include photon attenuation or scattering. 
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The conventional approach to the problem of collimator blurring assumes that the depth-dependent point-response 
function can be approximated by a spatially-invariant blurring function, so that inverse filtering can be applied to the 
projection data in each view independently 5 . Our approach involves simultaneous inverse filtering of the data in all views 
prior to conventional image reconstruction from projections, and is an extension of our earlier work 6 on properties of the 
Fourier decomposition of the sinogram. Ogawa et al? have recently presented a method for correction of depth-dependent 
collimator blurring that also involves Fourier decomposition of the sinogram, but their formulation of the problem and the 
coefficients of their filter function are different from ours. Fourier analysis of the sinogram has also been used by Hawkins 
et al % as the basis of an alternative reconstruction method for SPECT. Their formulation is directed toward solving the 
problem of photon attenuation (which is assumed to be uniform within a convex boundary), rather than the problem of 
depth-dependent collimator blurring which is the motivation for the preprocessing method in this paper. 


2. NOTATION AND PRELIMINARIES 


2.1 Coordinates and projections 

Our notation for object and detector coordinates in the 2-D plane is shown in Fig. 1. In object coordinates, (z, y) 
represents a point in the plane, and f(x,y) denotes the activity (emission density) at (z,y). The detector coordinate axes 
rotate with the collimator and detector about the origin as shown in Fig. 1, where the detector coordinates (s, t) are rotated 
by 9 from the fixed object coordinates (x,y). We denote by f*(s t t) the activity at (s,t) in rotated (detector) coordinates: 

f e (s,t) = f(scos9 -tsinO } ssind + tcos9) (1) 

A line in the plane parallel to the t axis is specified by its pair of parameters 0). We denote by p(s,9) the value of the 
line integral of activity along the line (s,9) : 

p(M) = J?(s,t)dt (2) 
The inversion process necessary to estimate / from p is known as image reconstruction from projections 9 ' 10 . 

2.2 Collimator response function 

The assumed collimator consists of a row of parallel holes, where the long axis of each hole is parallel to the t axis, 
and the row of holes is at a fixed distance t c from the center of rotation. In our 2-D model, these "holes" are actually 
1-D bins separated by partitions (septa). For a collimator hole centered on the line (s, 0), as shown in Fig. 2, the detected 
response to a point source of activity located at (s',t) in detector coordinates depends on t (the distance along the hole axis 
minus t c ) and on the distance s' - s from the point to the center-line of the hole. We denote this point-response function by 
ft, and we assume that the point response for a fixed t is symmetric about the center-line, so that 

h( s -*',*) = h(s' - M) = h(\s f - si t) (3) 
We denote by p c (s, 9) the value of the "fan integral" of activity for the collimator hole centered on the ray ($,9) : 

Pe (s,*) = J fis'^his-s'^ds' dt (4) 

In the following sections of this paper, we present a method to estimate p from p c so that conventional image reconstruction 
algorithms for p can be used for collimator projection data p c that has been preprocessed using this method. 


2.3 Point source geometry 

Consider a point source of activity at polar coordinates r p and </> pt as shown in Fig. 3. The (ar,y) coordinates of this 
point are {r p cos<j> p ,r p sin<j> p ) and the (s,<) coordinates are denoted by (s p (0),t p (9))* where 

Sp (9) = r p cos(<f> p - 9) and t p {0) = r p sin(</> p - 9) (5) 
The collimator response to this point source is 

Pc(s>0) = h(8-s p {9),l p (9)) (6) 
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2.4 Fourier analysis of point-response function 

The 1-D Fourier transform of h(s> t) with respect to its first variable is denoted by h(u, t) , where 

A(w,t) = [ + °° h(s,t)e^ s ds (7) 

J— oo 

For the shifted point-response function, we find 

h(«,t)e- ,W = h{s-s\t)e-^ s ds (8) 

*/ -oo 

2.5 Fourier analysis of the sinogram 

The sinogram is the 2-D map of p as a function of s and We make use of 2-D Fourier analysis of the sinogram, 
which is periodic in 9 (but not in s). The Fourier coefficients are denoted by p(w,n), where a; is a real number and n is 
an integer, and we express the Fourier pair as 

flu/,n) = — / / p(s,e)e-^*+ nd Usde (9) 

-I y + OO OO 

p( s ,0)=— Yl k",n)e«"" +n9) duj (10) 

^ J-OO _ 


3. FOURIER ANALYSIS FOR A POINT SOURCE 


3.1 Exact analysis 


The Fourier coefficients of the collimated projections are defined by 

Pc^.n) = -L pJs^e-^'WdsdB (11) 

£X Jo J-oo 

For the case of a point source at (r p cos$ p , r p sin<f> p ), we substitute h($ - s p (0),t p (0)) for p c (s } 0), then integrate over s 
to obtain 

i ft* 

p e (w, n) = — / /i(u;^ p (^))e- l ^ r ' co ^^-^ +ri ^ ^ (12 ) 

3.2 Asymptotic analysis 

We now examine the characteristics of the complex exponential factor in the integrand, and indicate how the principle 
of stationary phase 11 may be used to obtain approximations to the integral. Our purpose here is simply to outline the steps 
involved, without attempting to present a mathematically rigorous derivation. We introduce the notation 3>(0), where 

<b{0) = u>r p cos(<j> p - 0) + nO = ws p (Q) + nO (13) 

so that 

P c (u, n) = i- h{u, f p (0))e-'*W dO . . (14) 

Jo 

Figure 4 illustrates the behavior of $ and cos$ as functions of 0 for (w,n) = (-80,40). For w large, e - '*^ oscillates 
rapidly as a function of 0, except near angles (denoted by 6i and 0 2 ) at which b^/dO = 0. The relatively small portions 
of the range of integration centered around the angles 0i and 0 2 are the only ones that make significant contributions to the 


234 / SPIE Vol. 1092 Medical Imaging III: Image Processing (1989) 


# 


value of the integral, since the integral over the oscillatory parts of the function is almost zero. If h(u> y t p (0)) varies only 
slowly as a function of 0, compared with c" 1 '*^), then f** can be approximated by f£*+^ + f**** as follows: 


fc(«,n)« ^jf^Mw^pWJe-'^rftf + ^^(^pie^* (15) 
We now take h outside each of the integrals, since the intervals are small, and h is assumed to be slowly varying, so that 

We find $i and 0 2 from the condition for stationary phase (i.e., d$/dO = 0), where 


(16) 


do 


= vr p sin(<f> p - 0) + n = w« p (fl) + n (17) 


~- = 0 leads to sin(l - 0) = and tJO) = (18) 

00 wr p cj 

For \n\ < \ur P l there exist angles 0 lt 0 2 in the range [0,2ir] for which d^/80 = 0. We also require t p (0x) and < p (^) 
in the expression for ^ c (w,n), and the condition for stationary phase indicates directly that t p {9\) and t p (9 2 ) have the 
remarkably simple form 


t p (e 1 ) = t p (e 2 ) = -- t (19) 


so that 


= fc( w ,i p (fl 2 )) = A (w,- J) (20) 

and 

*(-■-)-*(-•-)' {sjC'r--"* 0 ** sjC*"*"} <2i) 

For the ideal case of the infinite-resolution collimator, p c — ► p and /i — ► 1, so that p(w, n) ^ { / + / } in the previous 
equation, and 

i> c (w,n) « P( w . n ) ( 22 ) 
This result shows that an estimate of p(u>,n) can be recovered, in principle, from p c (w, n) by 

j3(w,n) « [A(w,-n/w)] 1 p c (w,n) (23) 

where -n/w)] 1 denotes the inverse filter. 

The above Fourier analysis is based on the detected response to an arbitrary point source of activity. We find, however, 
that the final expressions for p c and p are independent of the coordinates of the point source. By the principle of superposition, 
and the linearity of the operations involved, these expressions are therefore applicable to the general case of a non-localized 
distribution of activity i 
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4. PROPERTIES OF THE FOURIER COEFFICIENTS 

4.1 Double- wedge characteristic 

The analysis of the Fourier coefficients of (he sinogram in the previous section shows that, for a point source of activity 
at polar coordinates r D and <j> p , the condition for stationary phase is satisfied at two angles within the range of integration, 
provided |n| < | W r p |. P Assuming r p non-negative, we have \n/u>\ < r p . Let R, be the largest value of r p f or any point in 
die activity distribution. Then R, is the largest value of \n/u>\ for which the stationary condition is satisfied for at least one 
point source in the activity distribution. For \n/u\ > R,, the stationary-phase condmon is satisfied only for points outside 
the activity distribution and, as a consequence, the Fourier coefficients are almost zero. Because of this property, the > Fourier 
coefficients form a double-wedge or "bow-tie" pattern in the the (w, n) plane, as shown in Fig. 5, which is described below. 
A more detailed discussion of the geometrical basis of this property of the Fourier coefficients is given in our earlier report , 
and we note in passing that the coefficients of other projection expansions involvmg angular harmonics are found to nave 
their own double-wedge properties, known as projection consistency conditions • . 

Figure 5 shows the results of experiments using a phantom consisting of two uniform disks of emission activity inside 
a unit circle, where the disks each have a radius of 0.1 and are centered at coordinates (0.8, 0) and (-0.8, 0). Figure 5a 
shows the sinogram of the simulated parallel projections p(s,0) from an ideal detector, the image reconstructed from this 
sinogram, and the map of the calculated Fourier coefficients j3( W , n). The double-wedge characteristic is evident Projection 
data were also simulated for the collimator case, where the point-iesponse function was calculated geomemcal y using die 
angle subtended at the point by the collimator aperture. Figure 5b shows the sinogram of the simulated collimator data, the 
image reconstructed from the collimator data, and the map of the Fourier coefficients of the collimator data. 

4.2 Connection with Bessel functions 

In the case of projections of a point source obtained by an infinite-resolution collimator, we find that die integral for the 
Fourier coefficients reduces to the integral representation for the Bessel function. Since the properties of Bessel functions 
have been the subject of extensive investigations 11 , many interesting properties of this integral, and hence properties of the 
Fourier coefficients of ideal data, can be deduced from them. 

The Bessel function of the first kind, order n, is denoted by J„, and ils integral representation is 

Jn( y) = J_ [ 2 * + " e i(ne-z« n ») de (24) 

where a is an arbitrary angle. For the ideal case of the infinite-resolution collimator, h — 1, so that the projections of a 
point source at polar coordinates r p and <j> P lead to 

'£( W)n ) = _L j'\-i[-r r co^ r -9)+ne\ M (25) 

2t J 0 

In terms of Bessel functions, £(w,n) may be expressed as 

^ w ,n) = J n ( W r p )e-«»^^ W 

For non-zero n, the Bessel function J n (z) has a different behavior for \z\ < \n\ than it has for | s | > |n|. * " > 1 • J«| 
for increasing values of the argument beginning from * = n, the function first increases to ^ p ^™^^«S5*K 
oscillates about the zero axis, similar to a decaying sinusoid, with a slowly decreasing amplitude of osculation asymptotically 
proportional to For decreasing values of beginning from z = n, the integral that represents the f ^sel funcuor idoes 
not include a point of stationary phase, and the function decreases rapidly and monotonically toward the zero axis, which 
it eventually intersects at z = 0. These characteristics of the Bessel function, together with 

|P(u,,n)|=|J„(wry)| ( 27 ) 

imply that the non-negligible values of p{w, n) are confined to a double-wedge shaped region in the ( w , n) W-J**™^ 
. n = ±u>r p mark the boundary of the double-wedge region. Considering a general object as a superposition of points, we still 
get the double-wedge shaped fi distribution as a result of the superposition of the transforms of the individual point objects. 
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The double- wedge shaped pattern in the (w, n) plane due to each point object has a transition region near the boundary 
lines n - ±u>r p . We now point out some relevant characteristics of the Bessel function in this transition region (see also 
Natterer 12 , pages 197-199). Let Z n denote the location of the first local maximum (this is also the global maximum) of 
J n (z). For n positive, and reasonably large (e.g., n > 10) the following large-order asymptotic approximations are useful 
(see Abramowitz and Stegun 13 , formulas 9.5.16, 9.5.20, 9.3.5 and 6.1.13): 

Z n -n + 0.809n* with J n (Z n ) ~ 0.674n"3 and J n (n) <- 0.447n~3 (28) 

For z < n, it is of interest to determine the value of z at which J n {z) is negligible compared to its maximum value. 
Numerically solving J n (z) = 0.017„(Z„) , we find that z ~ n - L35n3 . This value of z may be considered to mark the 
low- valued edge of the transition region. The width of the transition region is therefore dependent only weakly on n. 

4.3 Frequency-distance relation 

We first state in general terms a property of the Fourier coefficients that we call the "frequency-distance relation" 6 . We 
then show the implications of this property by means of an illustrative example. In section 3.2 we found that 

t,(*i) = « p (fl 2 ) = -» (29) 

This remarkable result is what we refer to as the frequency-distance relation for the Fourier coefficients. It means that as $ 
goes from 0 to 27r, the most significant contribution of an arbitrary point object to the value of the coefficient p(w y n) occurs 
at an angle value of 6 such that the point's distance t along the ray is -n/u. The importance of this relation lies in the 
possibility of modifying the coefficients within a narrow area around the line n — ~ut in the Fourier domain and obtaining 
the effect of modifying contributions to the projection data that have been emitted at a specific distance t from the detector. 
This is a new concept in projection imaging, initially presented in our previous report 6 , and it is not easy to present and to 
understand in the form of the general statement that we have just given. The discussion and specific examples that follow 
are intended to make clearer the meaning and implications of the frequency-distance relation for the Fourier coefficients. 

Since t in the detector coordinates represents a signed distance from the object location to the s axis, a negative value 
of t means that the point of interest is located in the half of the object which is closer to the detector. When t is of positive 
value the half of the object which is farther away from the detector is indicated. From the frequency-distance relation 
t p = — n/w, we find that in the (u>,r?) plane, the Fourier coefficients p in the two quadrants where n/cj > 0 correspond 
to contributions to the projection data from point sources that are nearer to the detector than the center of rotation (where 
t = 0). The other two quadrants corresponding to n/co < 0 represent contributions to the projection data from point sources 
that are farther from the detector than the center of rotation (Fig. 6). 

Assuming data are collected for a full 360° range of views, a particular point source in the object will have negative 
t (i.e., it will be nearer to the detector than is the center of rotation) for a subintervai of 180? within the full range of 
views. This subintervai of 180° is sufficient for accurate reconstruction of this point. Different point sources will have 
different subintervals, each of length 180°, located within the full range of views and corresponding to negative t in detector 
coordinates. The power of the Fourier analysis involving the frequency-distance relation is that it enables the full-range 
projection data to be decomposed into another data set where each point is represented by its own appropriate subintervai 
within the full range of views. 

A simple experiment illustrates this property. Wc see from the analysis above that if only the coefficients for which 
n/u) > 0 were selected, then the data used for reconstruction would involve only contributions from points near to the 
detector, where the resolution and sensitivity are best. The experiment involves the following steps. First, the emission 
projections p c (s,0) are obtained from the collimated detector model for a complete rotation in 0. Then, the 2-D discrete 
Fourier transform of p c (s,0) produces p c (uj,n). Now mask out the p c data from the second and fourth quadrants. Finally, 
the 2-D inverse Fourier transform of the selected coefficients in the first and third quadrants produces a modified sinogram 
p c (s,0). Reconstruction is done from these selected components of the projection data to get an improved image. 

Figure 7a shows the blurred sinogram from the simulated SPECT detector, its Fourier decomposition and the image of 
two blurred disks reconstructed from it From Fig. 7b we can see the improvement in the reconstruction after following the 
procedure mentioned above. In comparison, Fig. 7c shows what happens if we do the complementary procedure, selecting 
those quadrants that involve contributions to the projection data from point sources that are farther from the detector than is 
the center of rotation. The collimator resolution is poor for such point sources, as Fig. 7c shows. This experiment shows 
that heavily blurred data can be simply picked out and discarded from collimated projections by processing in the Fourier 
domain. The frequency-distance relation is more powerful than this, however, because it can be used to correct each blurred 
projection depth-by-depth along the projection ray without discarding any detected events, as we describe in the next section. 
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5. DEPTH-DEPENDENT SPATIAL FILTERING OF THE SINOGRAM 

5.1 Procedure for filtering of projection data 

The analysis in section 3.2 leads to the following Fourier method for correction of depth-dependent collimator blurring: 

1) Compute p c (u>, n) by the 2-D discrete Fourier transform of the data p c (s, 9) obtained for a complete rotation around 
the patient; _ 

2) Recover estimate of n) from £ e (w, n) by £(u>, n) « -n/w)] »). wherc lhe inverse filter is determined 
by the procedure described in section 5.2; 

3) Compute estimate of 0) by 2-D inverse Fourier transform. 

5.2 Inverse filter design 

The first step in computing the inverse filter is to specify the point-response function h($,t) of the collimator-detector 
combination. This combined response function has two components, namely, the geometrical response function of the 
collimator and the intrinsic point-response function of the detector. 

In experimental measurements of point-response functions, it has been found 14 that a Gaussian function in $ provides 
a good fit to the experimental data for each /. With a Gaussian representation for h(s,t), its 1-D Fourier transform 
can be found analytically for each t, so that h(u>, -n/w) can be evaluated at the required sample points in the (w, n) space. 

Once h(u, -n/w) is available, an inverse filter (which we denote by -n/u)] *) may be designed. The simple 
reciprocal of h(w,n) is a reasonable inverse filter only for those (w, n) where the signal-to-noise ratio is high and n)| 
is not too small. At high frequencies, the signal-to-noise ratio is low, and |h(u>, n)\ is small, so that the high-frequency 
characteristic of a useful inverse filter must deviate from the simple reciprocal in order to avoid excessive noise amplification. 
An appropriate high-frequency characteristic is also necessary in conventional data preprocessing in SPECT • 1 for the same 
reason. The Wiener filter and Metz filter have been used 16 " 19 to construct deblurring filters for SPECT projection data 
based on the assumption that the actual blurring can be approximated by a uniform blurring corresponding to the center 
of rotation. We are presently investigating the properties of depth-dependent inverse filters derived from Wiener, Metz, 
and other filter prototypes. 

5.3 Results 

Figure 8 shows the results of a preliminary simulation study to test the basic principles of the method. The phantom 
that was simulated consists of circular regions of various sizes, each having a uniform distribution of activity. Collimator 
datap c (s,0) was simulated for 135 samples in s, 360 samples in 0 (complete rotation) for a parallel-hole collimator having 
a hole depth/diameter ratio of 5. In this preliminary study, the data were noiseless (infinite statistics) with no attenuation or 
scattering. The geometrical point-response function of the collimator was used in the simulation of the data, and a Gaussian 
approximation was used for the purpose of calculating the inverse filter, as described in section 5.2. Although this simulation 
represents an idealized imaging situation, the result does indicate that the method has potential for deblurring projection 
data, and hence for improving the fidelity of the reconstructed image. 

6. CONCLUSIONS 

We conclude that a linear filtering approach based' on Fourier transforms has potential for the correction of depth- 
dependent collimator blurring in SPECT. In future work we plan to investigate the validity of the approximations introduced 
in the analysis, and to formulate the method for more general models of the data collection process. The design of a robust 
inverse filter for noisy data will be crucial for the successful application of the method to measured data in SPECT. 
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Fig. 1 Projection Geometry Fig.2 Geometry of coilimator response function for 

hole centered on line (s, 9) and point source at (a', t). 



Fig.3 A point object in both image Fig.4 Typical plots of $(d) and cos($(0)) for a fixed (w, n). 

(x, y) and detector (s, t) coordinates 
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Fig.6 Contribution 
of source points to 
the map of Fourier 
coefficients of the 
sinogram. 


a b 


Fig. 5 Comparison of the sinograms from an ideal 
detector (a) and a collimated detector (b), and 
corresponding reconstructed images from those 
sinograms, also the maps of Fourier coefficients of 
the sinograms. 
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a b c 


Fig J Demonstration of the simple simulation to modify SPECT sinogram and 
reconstruction, (a) is the collimated data, (b) shows the data modified by discarding 
Fourier coefficients which are corresponding to the part of projection data contributed 
by sources farther away from the detector than the center of rotation, (c) is the 
complementary part of the data of (b). It shows the blurring of collimator at large 
distances. 


242 / SPIE Vol. 1092 Medical Imaging III: Image Processing. (1989} 


7-B 


2-B 3-B 


Fig.8 Simulation of the data acquired from a collimated detector system. High sensitivity par- 
allel hole collimator with a bore depth to diameter ratio of 5:1 is simulated. Attenuation and septal pen- 
etration are ignored in the simulation. Deblurring filter is based on Gaussian shaped PSF assumption. 
Phantom is composed of circles with different sizes within an unit image area. Image size is 135 by 
135. Blurred sinogram is double Fourier decomposed and processed in Fourier domain, followed by an 
inverse Fourier transform to get a modified sinogram. The modified sinogram leads to a deblurred image 
reconstruction. 

1-A Fourier decomposition of the sinogram from an ideal detector system with distance 
. invariant sensitivity. 

1- B Reconstruction from the ideal sinogram by standard filtered back-projection algorithm. 

2- A Fourier decomposition of the blurred sinogram from the collimated detector system. 
2*B Reconstruction resulting from the blurred sinogram. 

3- A Fourier decomposition in 2-A processed by the deblurring filter. 

3-B Reconstruction resulting from the modified sinogram obtained by inverse transform of 
3-A. 
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